Differential expression
load("DataSummary/DESeq_Cohort1.RData")
options(stringsAsFactors = FALSE)
mTable = DESeq_Cohort1$rList[[3]]
mTable = data.frame(ENSG = rownames(mTable), mTable)
AnnT = DESeq_Cohort1$AnnT
head(mTable)
ENSG baseMean
ENSMUSG00000090025 Gm16088 ENSMUSG00000090025 Gm16088 2.013630e-01
ENSMUSG00000025900 Rp1 ENSMUSG00000025900 Rp1 7.937009e+00
ENSMUSG00000025902 Sox17 ENSMUSG00000025902 Sox17 2.084278e-02
ENSMUSG00000098104 Gm6085 ENSMUSG00000098104 Gm6085 2.618361e+00
ENSMUSG00000088000 Gm25493 ENSMUSG00000088000 Gm25493 5.372398e-02
ENSMUSG00000033845 Mrpl15 ENSMUSG00000033845 Mrpl15 1.402656e+03
log2FoldChange lfcSE stat
ENSMUSG00000090025 Gm16088 -0.028940390 0.09136916 -0.3167413
ENSMUSG00000025900 Rp1 0.106469263 0.14046002 0.7580040
ENSMUSG00000025902 Sox17 0.008723823 0.03637123 0.2398551
ENSMUSG00000098104 Gm6085 0.291415163 0.14016518 2.0790839
ENSMUSG00000088000 Gm25493 0.013212903 0.03637421 0.3632492
ENSMUSG00000033845 Mrpl15 0.056124421 0.02751339 2.0398950
pvalue padj
ENSMUSG00000090025 Gm16088 0.75143989 0.8628992
ENSMUSG00000025900 Rp1 0.44844856 0.6361525
ENSMUSG00000025902 Sox17 0.81044263 NA
ENSMUSG00000098104 Gm6085 0.03760964 0.1171989
ENSMUSG00000088000 Gm25493 0.71641875 NA
ENSMUSG00000033845 Mrpl15 0.04136079 0.1256043
cat(nrow(filter(mTable, padj < 0.01, log2FoldChange > 0)), " up-regulated genes.\n",
nrow(filter(mTable, padj < 0.01, log2FoldChange < 0)), " down-regulated genes.\n",
sep = "")
1592 up-regulated genes.
1652 down-regulated genes.
hTable = mTable[rownames(AnnT)[AnnT[, "Homolog_Human"] != ""], ]
cat(nrow(filter(hTable, padj < 0.01, log2FoldChange > 0)),
" up-regulated genes with human homologs.\n",
nrow(filter(hTable, padj < 0.01, log2FoldChange < 0)),
" down-regulated genes with human homologs.\n", sep = "")
1201 up-regulated genes with human homologs.
1361 down-regulated genes with human homologs.
Pathway enrichment
source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")
load("DataSummary/pwDB.RData")
pwDB <- pwDB[c(2:5,9,14:18)]
UP = unique(AnnT[filter(hTable, log2FoldChange > 0, padj < 0.01)$ENSG, "Homolog_Human"])
DN = unique(AnnT[filter(hTable, log2FoldChange < 0, padj < 0.01)$ENSG, "Homolog_Human"])
PG = unique(AnnT[hTable$ENSG, "Homolog_Human"])
pwDBList <- list()
for(i in 1:length(pwDB))
pwDBList[[i]] <- pwListFilter(pwDB[[i]], PG, minN = 15, maxN = 2000)
names(pwDBList) = names(pwDB)
write.table(as.matrix(pHyperListDB(UP, pwDBList, PG, FALSE)),
file = "out/rainbow_UP_PathwayEnrichment.txt",
row.names = F, col.names = T, sep = "\t")
write.table(as.matrix(pHyperListDB(DN, pwDBList, PG, FALSE)),
file = "out/rainbow_DN_PathwayEnrichment.txt",
row.names = F, col.names = T, sep = "\t")
cat(length(UP), " unique up-regulated genes with human homologs.\n",
length(DN), " unique down-regulated genes with human homologs\n", sep = "")
1200 unique up-regulated genes with human homologs.
1361 unique down-regulated genes with human homologs
uTable = read.table(file = "out/rainbow_UP_PathwayEnrichment.txt", row.names = NULL, header = TRUE, sep = "\t")
uTable = arrange(uTable, Pvalue)
uPathways = c("HALLMARK_MTORC1_SIGNALING",
"HALLMARK_OXIDATIVE_PHOSPHORYLATION",
"HALLMARK_MYC_TARGETS_V1",
"REACTOME_CHOLESTEROL_BIOSYNTHESIS",
"REACTOME_REGULATION_OF_MITOTIC_CELL_CYCLE",
"REACTOME_SIGNALING_BY_WNT",
"REACTOME_METABOLISM_OF_MRNA",
"REACTOME_ACTIVATION_OF_NF_KAPPAB_IN_B_CELLS")
uTable[uTable$PathwayName %in% uPathways, 1:7]
PathwayName DataBase
8 HALLMARK_MTORC1_SIGNALING MSigDB_h_all_v5
14 HALLMARK_OXIDATIVE_PHOSPHORYLATION MSigDB_h_all_v5
17 HALLMARK_MYC_TARGETS_V1 MSigDB_h_all_v5
24 REACTOME_CHOLESTEROL_BIOSYNTHESIS MSigDB_c2_cp_reactome
94 REACTOME_REGULATION_OF_MITOTIC_CELL_CYCLE MSigDB_c2_cp_reactome
97 REACTOME_ACTIVATION_OF_NF_KAPPAB_IN_B_CELLS MSigDB_c2_cp_reactome
111 REACTOME_SIGNALING_BY_WNT MSigDB_c2_cp_reactome
134 REACTOME_METABOLISM_OF_MRNA MSigDB_c2_cp_reactome
SignatureSize PathwaySize Overlap Pvalue Qvalue
8 1200 196 61 5.957365e-22 4.557384e-19
14 1200 190 54 1.196260e-17 5.229363e-15
17 1200 193 52 5.675148e-16 2.043053e-13
24 1200 20 15 9.558724e-15 2.437475e-12
94 1200 76 24 4.603195e-10 2.996974e-08
97 1200 61 21 6.883535e-10 4.343014e-08
111 1200 61 20 4.458776e-09 2.458352e-07
134 1200 188 38 2.843030e-08 1.293631e-06
dTable = read.table(file = "out/rainbow_DN_PathwayEnrichment.txt", row.names = NULL, header = TRUE, sep = "\t")
dTable = arrange(dTable, -Pvalue)
dPathways = c("HALLMARK_HEME_METABOLISM",
"HALLMARK_INTERFERON_GAMMA_RESPONSE",
"HALLMARK_INTERFERON_ALPHA_RESPONSE",
"HALLMARK_IL6_JAK_STAT3_SIGNALING",
"HALLMARK_KRAS_SIGNALING_UP")
dTable[dTable$PathwayName %in% dPathways, 1:7]
PathwayName DataBase SignatureSize
5792 HALLMARK_KRAS_SIGNALING_UP MSigDB_h_all_v5 1361
5852 HALLMARK_IL6_JAK_STAT3_SIGNALING MSigDB_h_all_v5 1361
5863 HALLMARK_INTERFERON_ALPHA_RESPONSE MSigDB_h_all_v5 1361
5971 HALLMARK_INTERFERON_GAMMA_RESPONSE MSigDB_h_all_v5 1361
6119 HALLMARK_HEME_METABOLISM MSigDB_h_all_v5 1361
PathwaySize Overlap Pvalue Qvalue
5792 188 29 1.695827e-03 3.154548e-02
5852 82 16 9.805450e-04 2.230831e-02
5863 81 16 8.477694e-04 2.010988e-02
5971 179 31 1.558219e-04 6.315430e-03
6119 180 71 3.090972e-29 9.458375e-26
uTable = uTable %>% mutate(group = "UP", LogFDR = -log10(Qvalue))
dTable = dTable %>% mutate(group = "DN", LogFDR = log10(Qvalue))
pwTable = rbind(uTable, dTable)
colorG = rep("NC", nrow(pwTable))
colorG[pwTable$LogFDR > -log10(0.05)] = "upSig"
colorG[pwTable$LogFDR < log10(0.05)] = "dnSig"
pwTable$colorG = colorG
pwTable = filter(pwTable, DataBase %in% c("MSigDB_h_all_v5", "MSigDB_c2_cp_reactome", "MSigDB_c2_cp_kegg", "MSigDB_c2_cp_biocarta"))
p <- ggplot(pwTable, aes(x = 1:nrow(pwTable), y = LogFDR, color = colorG))
p <- p + theme_bw() + labs(x = "Pathways", y = "-log10 FDR", title = "")
p <- p + theme(plot.title = element_text(hjust = 0.5)) + geom_point(size = 0.1, color = "grey")
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + geom_point(size = 0.2)
p <- p + scale_colour_manual(values = c("upSig" = "#dd1c77", "dnSig" = "#3182bd", "NC" = "grey"))
print(p)
Overlapping CRISPR hits
load("DataSummary/Focused_Depleted.RData")
upGS = unique(sapply(strsplit(filter(mTable, log2FoldChange > 0, padj < 0.01)$ENSG, " "), function(z) z[2]))
dnGS = unique(sapply(strsplit(filter(mTable, log2FoldChange < 0, padj < 0.01)$ENSG, " "), function(z) z[2]))
pgGS = unique(sapply(strsplit(mTable$ENSG, " "), function(z) z[2]))
BG = intersect(pgGS, PG)
UP = intersect(upGS, BG)
DN = intersect(dnGS, BG)
### MA cell line
res <- matrix(NA, length(maList), 6)
for(i in 1:length(maList)){
Depleted = intersect(maList[[i]], BG)
U1N = length(intersect(UP, Depleted))
U0N = length(setdiff(UP, Depleted))
D1N = length(intersect(DN, Depleted))
D0N = length(setdiff(DN, Depleted))
ABCD = matrix(c(U1N, U0N, D1N, D0N), nrow = 2, ncol = 2, byrow = TRUE)
dimnames(ABCD) = list(c("UP", "DN"), c("Depleted", "NotDepleted"))
res[i, 1:4] = c(U1N, U0N, D1N, D0N)
res[i, 6] = fisher.test(ABCD)$p.value
res[i, 5] = fisher.test(ABCD)$estimate
}
colnames(res) = c("UP_Depleted", "UP_NotDepleted", "DN_Depleted", "DN_NotDepleted", "OR", "Pvalue")
rownames(res) = names(maList)
knitr::kable(data.frame(res))
| Post_vs_plasmid |
102 |
23 |
39 |
22 |
2.4882466 |
0.0107318 |
| Post_vs_Pre |
92 |
33 |
29 |
32 |
3.0559849 |
0.0005989 |
| Post_vs_vitro |
15 |
110 |
13 |
48 |
0.5054908 |
0.1255157 |
| Pre_vs_plasmid |
100 |
25 |
31 |
30 |
3.8390473 |
0.0000710 |
| vitro_vs_plasmid |
98 |
27 |
31 |
30 |
3.4860765 |
0.0001872 |
| vitro_vs_Pre |
93 |
32 |
29 |
32 |
3.1847263 |
0.0004950 |
### HM cell line
res <- matrix(NA, length(hmList), 6)
for(i in 1:length(hmList)){
Depleted = intersect(hmList[[i]], BG)
U1N = length(intersect(UP, Depleted))
U0N = length(setdiff(UP, Depleted))
D1N = length(intersect(DN, Depleted))
D0N = length(setdiff(DN, Depleted))
ABCD = matrix(c(U1N, U0N, D1N, D0N), nrow = 2, ncol = 2, byrow = TRUE)
dimnames(ABCD) = list(c("UP", "DN"), c("Depleted", "NotDepleted"))
res[i, 1:4] = c(U1N, U0N, D1N, D0N)
res[i, 6] = fisher.test(ABCD)$p.value
res[i, 5] = fisher.test(ABCD)$estimate
}
colnames(res) = c("UP_Depleted", "UP_NotDepleted", "DN_Depleted", "DN_NotDepleted", "OR", "Pvalue")
rownames(res) = names(hmList)
knitr::kable(data.frame(res))
| Post_vs_plasmid |
102 |
23 |
36 |
25 |
3.0589979 |
0.0013226 |
| Post_vs_Pre |
99 |
26 |
38 |
23 |
2.2935265 |
0.0205806 |
| Post_vs_vitro |
12 |
113 |
8 |
53 |
0.7049096 |
0.4607979 |
| Pre_vs_plasmid |
95 |
30 |
25 |
36 |
4.5176767 |
0.0000042 |
| vitro_vs_plasmid |
95 |
30 |
37 |
24 |
2.0456464 |
0.0388821 |
| vitro_vs_Pre |
94 |
31 |
31 |
30 |
2.9158955 |
0.0014460 |
sessionInfo
R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.6 (Maipo)
Matrix products: default
BLAS/LAPACK: /sibcb2/bioinformatics/software/Miniconda3/lib/libopenblasp-r0.3.15.so
locale:
[1] C
attached base packages:
[1] parallel grid stats graphics grDevices utils
[7] datasets methods base
other attached packages:
[1] viper_1.26.0 survcomp_1.42.0 prodlim_2019.11.13
[4] gridExtra_2.3 survHD_0.99.1 survC1_1.0-3
[7] Hmisc_4.5-0 Formula_1.2-4 lattice_0.20-44
[10] penalized_0.9-52 survival_3.2-11 Biobase_2.52.0
[13] BiocGenerics_0.38.0 stringr_1.5.0 dplyr_1.1.2
[16] scales_1.2.0 RColorBrewer_1.1-3 pheatmap_1.0.12
[19] VennDiagram_1.6.20 futile.logger_1.4.3 ggrepel_0.9.1
[22] ggplot2_3.4.2 rmarkdown_2.20
loaded via a namespace (and not attached):
[1] segmented_1.3-4 survivalROC_1.0.3 nlme_3.1-152
[4] tools_4.1.0 backports_1.2.1 bslib_0.2.5.1
[7] utf8_1.2.2 R6_2.5.1 KernSmooth_2.23-20
[10] rpart_4.1-15 rmeta_3.0 DBI_1.1.3
[13] mgcv_1.8-36 colorspace_2.0-3 nnet_7.3-16
[16] withr_2.5.0 tidyselect_1.2.1 downlit_0.4.2
[19] compiler_4.1.0 cli_3.6.1 formatR_1.11
[22] htmlTable_2.2.1 labeling_0.4.2 bookdown_0.33
[25] sass_0.4.0 checkmate_2.0.0 proxy_0.4-27
[28] digest_0.6.29 mixtools_1.2.0 foreign_0.8-81
[31] base64enc_0.1-3 jpeg_0.1-9 pkgconfig_2.0.3
[34] htmltools_0.5.2 parallelly_1.32.0 fastmap_1.1.0
[37] highr_0.9 htmlwidgets_1.5.4 rlang_1.1.1
[40] rstudioapi_0.15.0 SuppDists_1.1-9.7 jquerylib_0.1.4
[43] farver_2.1.1 generics_0.1.3 jsonlite_1.8.7
[46] distill_1.5 magrittr_2.0.3 Matrix_1.3-4
[49] Rcpp_1.0.9 munsell_0.5.0 fansi_1.0.3
[52] lifecycle_1.0.3 stringi_1.7.8 yaml_2.3.5
[55] MASS_7.3-54 listenv_0.8.0 splines_4.1.0
[58] knitr_1.42 pillar_1.9.0 future.apply_1.9.0
[61] codetools_0.2-18 futile.options_1.0.1 glue_1.6.2
[64] evaluate_0.20 latticeExtra_0.6-29 lambda.r_1.2.4
[67] data.table_1.14.2 png_0.1-7 vctrs_0.6.3
[70] bootstrap_2019.6 gtable_0.3.0 kernlab_0.9-29
[73] future_1.26.1 cachem_1.0.6 xfun_0.37
[76] e1071_1.7-11 class_7.3-19 tibble_3.2.1
[79] memoise_2.0.1 cluster_2.1.2 lava_1.6.10
[82] globals_0.15.1